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Abstract 

The hard sphere model is known to show a liquid-solid phase transition, with the 
solid expected to be either face centered cubic or hexagonal close packed. The 
differences in free energy of the two structures is very small and various attempts 
have been made to determine which structure is the more stable. We contrast the 
different approaches and extend one. 
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I. Introduction. 

Computer simulations of the hard sphere model in the late 1950's (see [B] for a 
review) showed a liquid/solid phase transition at packing fractions in the range 0.49 
to 0.54. Ever since that time there has been significant interest in determining the 
internal structure of that solid phase, since the model is one of the simplest possible 
showing both liquid and solid phases: the only interaction is that the spheres may 
not overlap. At least at high density, the solid is expected to have one of the 
structures associated with densest packing (packing fraction tt/vTS ~ 0.74), in 
particular the face centered cubic (fee) and hexagonal close packed (hep) structures, 
which commonly appear in materials. In 1968 Stillinger and coworkers [RSYS] used 
series expansions to try to answer the much simpler question of whether fee or hep is 
more stable. They computed the first 3 terms in a series for the entropy per sphere 
for both fee and hep and concluded that hep had higher entropy per sphere and 
was therefore more stable; no error bars were given. We will extend this analysis 
to 2 more terms. Others used molecular dynamics or Monte Carlo simulations 
[HR,FL] to compare the fee and hep entropies, with inconclusive results. In 1997, 
using more advanced computer technology. Woodcock [W1,W2] did a molecular 
dynamics simulation and concluded that fee is more stable, with error bars to back 
up this conclusion. The error bars and the magnitude of the entropy difference 
were called into question, but further simulations [BFMH,BWA,MH] agreed with 
Woodcock's qualitative claim that fee is more stable than hep. Finally, a geometric 
approach was recently applied to the first nontrivial term in Stillinger's expansion 
[RS]. This paper concerns the effort to extend the computation of the Stillinger 
expansion so as to determine the geometric mechanism behind stability. 

II. Analysis 

We consider the hard sphere model of classical statistical mechanics in the 
canonical ensemble, consisting of the uniform distribution of configurations of non- 
overlapping spheres at fixed packing fraction. (The momentum variables can be 
integrated out, as they are decoupled from physical space variables in this model.) 
We assume, based on old simulations [B], that at high density (packing fraction 
d ~ tt/vTS ~ 0.74) the model is in a solid phase which is, moreover, either an 
fee or hep crystal. More specifically, the assumption is that the distribution gives 
overwhelming weight to configurations of spheres which are perturbations of either 
a perfect fee or hep structure at the appropriate density. We begin with an analysis 
of the meaning of "perturbation" and of error estimates for this problem. 

One way to make definite the meaning of perturbations of the perfect fee struc- 
ture for d tt/vTS is as follows. Start with a finite block of spheres frozen in a 
perfect fee arrangement at density below close packing (obtained by relative choice 
of sphere radius and lattice spacing), and with a subset Q of \Q\ = N lattice 
sites for which we free up the associated spheres. As in [RSYS] we note that, as 
d 7r/\/l8, the restrictions in position of any sphere due to its neighbors can be 
approximated by neglecting the curvature of the spheres, obtaining linear restric- 
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tions. Specifically, suppose our lattice spacing is 1, Lj represents the coordinates of 

the J*'' lattice site, and Uj represents the relative coordinates of the center of the 
sphere of diameter 1 — e associated with this site. For nearest neighbor sites j and 
k, the hard core restriction is \\Lj + uj — — Uk\\ > 1 — e and the linearization 
is {Lj + Uj — Lk — Uk) ■ {Lj — Lfc) > 1 — e, or {uj — uu) • {Lj — Lk) > — e. So the 
entropy density s of this class of configurations at high density is precisely lia.{V)/N, 
where V is the volume of the convex polyhedron in R*^^ whose faces are the linear 
constraints on the positions of the N moveable spheres due to the frozen spheres 
and each other. 

Now that we have a class of configurations associated with each of fee and hep 
at high density, we consider the corresponding entropy densities s/cc and Shcp- As 
shown in [RSYS], Sj{d) diverges as d ^ tt/vTS and is of the form ln{Ad) + Cj + 
0(A(i), where Ad = n/^/TS — d and the Cj are constants. In order to determine 
whether fee or hep is more stable one must estimate the difference AC — Chcp—C fee- 
All works we are considering claim that |AC| fa 10~^. Using V for phase space 
volume we have by definition Vj = e^^^ , so obtaining the desired accuracy in Sj to 
order 10""^ for each case requires a bound on the relative error for the volumes of 
AVj/Vj <e^/^^' -1. 

A different way to make sense of perturbations of the fee and hep structures is 
to consider deviations of these structures that are periodic with some fixed period; 
if the density is high enough the two sets of configurations are disconnected and 
thus the volumes V/cc and Vhcp are well defined. Both molecular dynamics [W1,W2] 
and Monte Carlo simulations [BFMH,BWA,MH] used this setup, although mostly 
at densities near melting, or equivalently, for spheres of radius about 0.9, where 
the two sets of configurations are actually connected to one another (for periods 6 
or larger in each direction). At these densities, a simulation should theoretically 
sweep the entire phase space, regardless of where it is started. However, on the 
time scale available in practice the process is found to be trapped in regions whose 
configurations can be associated with a single lattice structure, and it is the volumes 
of these regions that have been used to conclude that fee is the more stable structure. 
Although the (extremely large) fcc-hcp mixing time may not be a problem in these 
simulations, the success of the method docs depend on the absence of relevant time 
scales beyond the ones that have been observed so far. Moreover, such simulations 
do not give any intuition into the mechanism, almost certainly geometric, which 
distinguishes these close packed structures. 

We now consider the method of Stillinger et al, which does not (inherently) 
rely on simulation and does in principle allow for a geometric interpretation. The 
first nontrivial term in the expansion below has been analyzed in [RS]. Consider 
again our finite block of spheres frozen in a perfect crystal arrangement, fee or hep, 
at density below close packing, and a subset Q of them which we free up. The 
entropy Sq of the moveable spheres associated with (or labelled by) the sites Q can 
be written as an exact sum 

PQQ PCQ:|P| = 1 PCQ:|P|=2 
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which defines the correction terms Cp. For each n > there are, up to symmetries, 
a fixed number of contributing "polymers" P of \P\ = n spheres; one can compute 
these and then each "n-body correction" Sn,Q = ^pcQ-\p\=n^P^ that 

'^Q \ ^ 

n>l 

As |Q| — > oo the number of occurences per site of the polymer P approaches a fixed 
frequency /p, and s„ q approaches a limit = J2\p\=n fp ^P- '^^^ entropy per 
sphere of the hard sphere solid is then s — ^„>i Sn- 

We specialize to densities very close to the maximum density, that is, phase 
space volumes (and hence entropies) are computed only to leading order in the 
difference e between the lattice spacing 1 and the sphere diameter 1 — e. As noted 
above exp(5'p) is the volume of a polyhedron in 3|P| dimensions, and the dimensions 
of this polyhedron all scale as e. In particular for \P\ — 1, si — Cp is simply the 
logarithm of the scaled volume of a Voronoi cell for a sphere, and this value is the 
same for fee as for hep. Note that for |P| > 1 the (divergent) ln(e) terms in Cp 
cancel, so s„ is well-defined in the e — > limit for n > 1. 

It is natural to simply truncate the sum s„ and compute each Sn by com- 
puting the appropriate volumes. See Table 1, below. The success of this method 
depends on the rapid decay of the individual terms Sn with increasing n, so that the 
sum of uncomputed terms are convincingly negligible. Note that for both fee and 
hep, and for all computed levels, is in absolute value considerably smaller than 
Sn- Recalling that we want to estimate As to within an error of 10~^, and assuming 
the two series remain geometrically decreasing by a factor of roughly 3, one needs 
to compute all terms up to level n= 7. Unfortunately this seems unattainable with 
current computers as we discuss in the next section. 
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-0.1647410 
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-0.0010935 
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-0.0517903 


0.8820808 


-0.0518249 


0.8831398 


-0.0010589 


4 


-0.0274325 


0.8546482 


-0.0261820 


0.8569577 


-0.0023094 


5 


0.0060605 


0.8607088 


0.0026023 


0.8595601 


0.0011486 



Table 1. Partial entropies for fee and hep structures. 

In addition to computing and adding terms in order, we can examine the contri- 
butions of individual polymers. Figure 1, below, shows the distribution of entropy 
corrections Cp that contribute to S5, after convolution by a Gaussian with variance 
lO"^'^. Note that the contribution of individual polymers should not be viewed as 
independent random variables. The sum of the ~ 500 terms contributing to S5 for 
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hep is not of order V500 larger than a typical term: rather, it is the same size as a 
typical term. Likewise for S3, S4, and for fee. Hopefully, by studying the distribu- 
tion of terms contributing to S3, S4 and S5 we will be able to devise resummation 
schemes with accelerated convergence. 
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Figure 1. Distribution of entropy corrections Cp contributing to S5. 



Another important goal is to understand the geometric mechanisms underlying 
the difference in entropy between fee and hep. The qualitative differences between 
the terms contributing to S2 were considered in [RS]. Although some general pat- 
terns have emerged from the data presented here (e.g., compact polymers tend to 
contribute negatively, while extended polymers tend to contribute positively), a 
true geometric understanding still eludes us. 



III. Computation 

The computations by Stillinger et al of (s2 and) S3 to 5-digit accuracy was 
remarkable for the time. We have taken advantage of current computers to compute 
Sn up to n = 5. To obtain S5 we needed to compute among other things the volumes 
of 591 polyhedra in 15 dimensions (the phase space for 5 spheres). The basic 
algorithm is simple, based on successive "triangularizations" : starting with d = 15, 
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choose a point in the (i- dimensional polyhedron to serve as the common apex of 
pyramids with the faces of the polyhedron as bases, and then add up the pyramid 
volumes AB/d, where A is the altitude and B the volume of the base. Given 
that the base is again a polyhedron, now in dimension d — 1, the process can be 
iterated. The main drawback of this algorithm is that the same lower dimensional 
volumes are computed repeatedly in different branches of the rcciusion. In the 
case of our 15-dimensional polyhedra, which have on average 55 faces, this leads 
to a highly impractical number of computations. The standard way of dealing 
with this problem is to store and to re- use the results from previous computations. 
Programs that implement this strategy are readily available, but for reasons of 
efficiency (the program that was recommended to us was too slow by several orders 
of magnitude) we wrote the necessary code ourselves; it is available at xx. Volumes 
in 15 dimensions each took between 7 and 160 hours of computation time, with 
an average of 29 hours, on a 64-bit processor running at 2.2 GHZ, using 2GB of 
memory. The number of re-used lower dimensional volumes in each of these cases 
is of the order of 10^°. Graphs of the entropy corrections fpCp associated with the 
volumes P are shown in Fig 1. 

We do not expect such an exact computation of sq to be practical before sub- 
stantial progress in computer hardware. 
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